/*
 * src/mpreal_friends.cpp
 *
 * This work was supported by the Director, Office of Science, Division
 * of Mathematical, Information, and Computational Sciences of the
 * U.S. Department of Energy under contract number DE-AC03-76SF00098.
 *
 * Copyright (c) 2004
 *
 */
#include <cfloat>
#include <arprec/mp_real.h>
#include "small_inline.h"

using std::cerr;
using std::endl;

/*
 * This Evaluates the gamma function,
 * using an algorithm of R. W. Potter
 * August 10th, 2005
 */
mp_real_temp gamma(const mp_real& t) {
    int i, j, ndp, neps, nt, nwords;
    double alpha, d2;
    static const double con1 = 1.151292547;
    mp_real eps, sum1, sum2, t1, t2, t3, t4, tn, z;

    nwords = mp::prec_words;

    // Is the following trucation ok ?
    neps = static_cast<int>((-nwords - 1) * 14.44943979);
    ndp = static_cast<int>((nwords - 1) * 14.44943979);

    eps = pow(mp_real(2.0), (-48*nwords - 48));

    // Handle special arguments.
    if (abs (t) > 1.0e8) {
	cerr <<"gamma: agrument is too large" << endl;
	return mp_real(-1.0).toTempAndDestroy();
    } else if (t == anint (t) && t <= 0.0) {
	cerr <<"gamma: invalid negative argument" << endl;
	z = mp_real(0.0);
	return z.toTempAndDestroy();
    }

    if (t > 0.0) {
	nt = static_cast<int>(dble(t));
	if (t == nt) nt = nt - 1;
	t1 = mp_real(1.0);

	for (i = 1; i <= nt; ++i)
	    t1 = t1 * (t - i);

	tn = t - nt;
	if (t == aint (t)) {
	    z = t1;
	    return z.toTempAndDestroy();
	}
    } else {
	nt = static_cast<int>(dble(1.0 - t));
	t1 = mp_real(1.0);
	
	for (i = 0; i <= nt - 1; ++i)
	    t1 = t1 / (t + i);
	
	tn = t + nt;
    }
    
    // Calculate alpha, then take the next highest integer value, so that
    // d2 = 0.25 * alpha^2 can be calculated exactly in double precision.
    
    alpha = aint (con1 * ndp + 1.0);
    t2 = tn;
    d2 = 0.25 * alpha * alpha;
    t3 = 1.0 / t2;
    sum1 = t3;

    // Evaluate the series with t, terminating when t3 < sum1 * epsilon.
    
    for (j = 1; j <= 1000000000; ++j) {
	t3 *= d2 / (j * (t2 + j));
	sum1 += t3;
	if (abs (t3) < abs (sum1) * eps) break;
    }
    if ( j == 1000000001 ) {
	cerr << "gamma: loop overflow 1" << endl;
	sum1 = 0.0;
    }
    sum1 = t2 * pow(0.5 * alpha, t2) * sum1;

    t2 = - tn;
    t3 = 1.0 / t2;
    sum2 = t3;

    // Evaluate the same series with -t, terminating when t3 < sum2 * epsilon.

    for (j = 1; j <= 1000000000; ++j) {
	t3 *= d2 / (j * (t2 + j));
	sum2 += t3;
	if (abs (t3) < abs (sum2) * eps) break;
    }
    if ( j == 1000000001 ) {
	cerr << "gamma: loop overflow " << endl;
	sum2 = 0.0;
    }
    sum2 = t2 * pow(0.5 * alpha, t2) * sum2;

    // Conclude with this square root expression.

    z = t1 * sqrt(mp_real::_pi * sum1 / (tn * sin(mp_real::_pi * tn) * sum2));
    return z.toTempAndDestroy();
}

